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Abstract 

We present a method that formally calculates exact frequency shifts of an electro- 
magnetic field for arbitrary changes in the refractive index. The possible refractive 
index changes include both anisotropic changes and boundary shifts. Degenerate 
eigenmode frequencies pose no problems in the presented method. The approach 
relies on operator algebra to derive an equation for the frequency shifts, which 
eventually turn out in a simple and physically sound form. Numerically the equa- 
tions are well-behaved, easy implement able, and can be solved very fast. Like in 
perturbation theory a reference system is first considered, which then subsequently 
is used to solve another related, but different system. For our method precision is 
only limited by the reference system basis functions and the error induced in fre- 
quency is of second order for first-order basis set error. As an example we apply our 
method to the problem of variations in the air-hole diameter in a photonic crystal 
fiber. 

Key words: 

PACS: 41.20.-q, 41.20.Jb, 42.25.-p, 42.81.-i 



1 Introduction 



Scale invariance of Maxwell's equations provides unique and simple relations 
between the electromagnetic properties of dielectric structures which only dif- 
fer by a simple spatially uniform scaling of the dielectric function [1]. The 
frequency shifts of electromagnetic modes due to spatially non-uniform low- 
index-contrast changes have typically successfully been addressed with pertur- 
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bation theory (see e.g. Ref. [2]) being quite similar to the one used widespread 
in quantum mechanics (see e.g. Ref. [3]). However, the development of fields 
such as photonic band-gap materials [1,4] and surface-plasmon subwavelength 
optics [5,6] calls for an understanding of high- index-contrast dielectric varia- 
tions beyond standard perturbation theory. There has recently been an effort 
toward a perturbative description of high-index-contrast problems with shift- 
ing material boundaries [7,8,9,10]. In Ref. [7] the problem is addressed by 
approaching the strong dielectric contrast as a small perturbation in the ma- 
terial boundary. This different approach to standard perturbation theory uses 
additional ideas from quantum mechanics, such as multiplying the perturb- 
ing dielectric function by a coupling strength variable. The related problem 
of scattering of electromagnetic waves by high-index-contrast dielectric struc- 
tures has been addressed by essentially non-perturbative approaches relying on 
iterative solutions of the Dyson equation and the Lippmann-Schwinger equa- 
tion for the electromagnetic Green's function and the electromagnetic fields, 
respectively [11,12,13]. 

Here, we describe a non-perturbative approach to electromagnetic eigenmodes 
that formally calculates exact frequency shifts of an electromagnetic field for 
arbitrary changes in the refractive index. As an example we apply our theory 
to photonic crystal fibers [14] and show how the approach can be used to study 
the effect of variations in the air-hole diameter. 

The paper is organized as follows: In Sec. 2 we present the formalism and in 
Sees. 3 and 3 we present the variational method in electric and magnetic field. 
In Sec. 5 we apply our theory to photonic crystal fibers. Finally, in Sec. 6 
discussions and conclusions are given. 



2 General formalism 



2.1 The temporal harmonic problem 



In the following we consider temporal harmonic solutions to Maxwell's equa- 
tions. The electrical field and magnetic field may then be written as 



E(r,t) = E,(r)e--*, (1) 
H(r,i) = H,(r)e--*, (2) 

with the subscript uj indicating the frequency dependence. The electrical field 
is a solution to a vectorial wave equation, which has the form of a generalized 
eigenvalue problem [1,15], 
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V X V X E,(r) = eo£(r) i^-j E^(r), (3) 

where e is the 3x3 dielectric tensor. The corresponding eigenvalue problem 
in the magnetic field has the form [1,15] 



V X eo^e-^(r)V x H,(r) = (^)'H,(r). (4) 

We consider lossless dielectrics so that e is real and symmetric positive semi- 
definite, i.e. it has non-negative eigenvalues. In the following we will for sim- 
plicity put c = eo = 1 and omit the subscripts on E^^ and H^. 



2. 2 Operator formalism 



In analogy with quantum mechanics we will in this work take advantage of 
a function space known as a Hilbert space, which consists of a set of basis 
functions and an inner product. In the following we use the Dirac bra-ket 
notation and define the inner product as 



A^)^ J dv Alir) ■ A^{v) 



(5) 



where f is the Hermitian conjugate (transpose and complex conjugate) of the 
field vector. If the basis of the Hilbert space is complete then any function in 
the Hilbert space can be written as a linear combination of basis functions, 
thus we can rewrite any operator O as 



where O has matrix elements 



(6) 



= /dr AL(r)0(r)A„(r). (7) 

Typically, the basis is not complete, but in many practical applications this is 
not imperative. Note, that Vx and i are no longer a differential operator and 
a tensor, respectively, but linear operators in the Hilbert space. 
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2.3 The reference system 



We will in the following consider a system that is similar to the systems we are 
trying to solve. We refer to it as the 'reference' system, and it is described by a 
dielectric system with dielectric function Eq, and its modes all have frequency 
uiq, but possibly with different wave vectors k . The corresponding operators 
will be assigned a subscript '0'. The eigenmodes are determined by either 
Eq. (3), or Eq. (4). 



3 Variational approach in E 



From Eq. (3) it follows that the frequency of the field \Ej can be determined 

by 





V X V X 






e 





(8) 



and since both V x Vx and i are Hermitian this forms the typical starting 
point for variational calculations. Here, we will show a different route which 
is beneficial for electromagnetic modes with the same frequency, but possibly 
with different wave vectors k. 

We may choose the eigenmode basis of the reference system, <En> , where 

I- J n=l 

eigenmodes are orthogonal with the norm 



En) = 5mn. (9) 

Here, Eq is the dielectric function operator and 5 is the Kronecker delta func- 
tion. 

In the spirit of perturbation theory we now solve a different but assumed 
similar system using the reference eigenmodes and frequencies. In quantum 
mechanics such a system is for historical reasons often referred to as the in- 
teracting system though the approach of course also can be applied to one- 
particle problems. Here we will decline from using the term interacting system 
and rather refer to it as the new system, since we seek to avoid linking our 
method to approximate perturbation theory. For the new system we denote 
the eigenfrequencies by uj and the new system also obeys the eigenvalue equa- 
tion, Eq. (3), and is described by the dielectric function, s(r). It is important 
to note that the frequencies of the new system do not need to be identical. In 
the new basis (eigenstate of reference system n) we define 
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V X V X 



En). 



(10) 



The double curl operation acting on a reference system eigenmode can be 
replaced by Eou^q and the equation rewritten as 



io(2}Q\En) = ifl En) ^ EqQ 



En) — ^0^ ^^O'-^o 



Er 



(11) 



where we have replaced the operator uj^ by the corresponding eigenvalue uj^. 
Next, we define 



T = £q£ — £q 



2 ' 



(12) 



so that our problem, Eq. (3), can be rewritten as 

E) = £qu'^\E] 



(13) 



Obviously, UqT is a Hermitian operator on our Hilbert space since (^ioi ^£o)^ — 
£o£~^£o. Similarly, igcu'^ is also Hermitian since ig is Hermitian. 

3.1 The variational problem 



We know expand the new system eigenmode(s) in the basis of reference system 
eigenmodes. 



E)^^Cn\E, 



(14) 



and since we have Hermitian operators on both the left-hand and right-hand 
sides of Eq. (13) we get the variational problem 



E 



cUq mm ■ 





f 











lin^ I 



'■^0 ^^^^2^CmCn{E„^ 



Er 



(15) 



where c„'s are normalized so that J2n l^nP = 1- It is of course interesting to 
know how much error in frequency an error in the state vector induces. For 
the numerical implementation this is particular important since here the basis 
set is rarely complete. Due to the variational principle, a first order error in 
the eigenmode E(r) will induce only a second error in frequency cu. 
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3.2 Finding uj 



We now want to find the eigenfrequencies of the new system, which we denote 
ijjm, m — 1,2, . . . ,N . First we find the matrix elements of the operator Eq. (11) 
as 



T 

J- rr 



E„ 



drEt^(r)£o(r)s-i(r)£o(r)E„(r) 
= / drDt^(r)£-^(r)D„(r), 



(16) 



which may be expressed both in terms of E and D fields. The variational 
problem Eq. (15) can be solved by finding the eigenvalues and eigenvectors of 
the matrix ui'^T. Formally we thus have 



Ci^^C^u^C (17) 
and the 'new' eigenfrequencies are given by 



u}^\/c(l^C'i ^CCIC\ (18) 

where a) is a diagonal cigenfrequency matrix and C* is a unitary transformation 
matrix, (7(7^ — 1. The 'new' system eigenstate vectors now become 



N 



m=l 



Being interested in the frequency shifts rather than the new frequencies, one 
subtracts the reference frequencies from the new frequencies defined by the 
eigenvalues of the operator given by Eq. (17). This gives 



ALJrn='^m- i^O- (20) 

It should be emphasized that in general mode m in the reference system may 
not be similar to the mode m in the new system. This is due to the fact 
that changes in the refractive indices may favor one mode more than another, 
causing different effective refractive indices of the modes to converge and then 
diverge again with the difference now having opposite sign of before. 
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3.3 Relation to first- order perturbation theory 



Consider a small variation in the dielectric function, 5e — s — Eq. This variation 
is assumed small enough to be neglected beyond first order. If the operators 

5e and e~'^ commute, = 0, (a sufficient, but not necessary, condition 

for this commutator to be zero is the two dielectric tensors being isotropic) 
the elements of the square frequency operator ui^T are 



Si 



(21) 



Expanding (l = ujq\I T"^ in 5e = Q gives 



K) + C»(fe') (22) 

where the new frequencies are the eigenvalues of O. This is the well-known 
result of first-order (degenerate) electromagnetic perturbation theory (see e.g. 
Refs. [2,7]). 



O — A _ Zll W 



Si 



4 Variational approach in H 



The H-ficlds satisfy eigenvalue equation Eq. (4) and all modes of the reference 
system again share the same frequency ujq, and satisfy the orthonormality 
Hnj — Smn- Next, we consider a presumably similar system 
dielectric function, and expand the modes of the 'new' system 



condition {^Hm 
with a different 



in the basis functions of the reference system: 



G 



V X £-V X H^) = Vx fi'W X 



Hn), 



(23) 



where wc have introduced the operator / = i~^io. Taking the dielectric func- 
tions to be anisotropic f(r) = l3x3£~^(r)£o(i") and using vector analysis we 
may rewrite the expression as 



V X /£o V X Hn) = V/ X £o V X Hn) + /V X V X 

^-iujVfx En) + fu^Hn), 



Hr 



(24) 
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and the eigenvalue equation (23) thus reduces to the operator equation where 
© has matrix elements 



Qrnn = -«^^o(^^m V/ X En) + U;^{Hm f Hn). 



(25) 



Since is a Hermitian operator we may apply the variational principle 



H 





e 









mm 



y'^i mn 



ioJolHm V/ X 



Er 



+ ^0 



f 



Hr, 



(26) 



Thus the variational approach in H gives a second order error in frequency 
for a first order error in the basis function. Note, that although the E field 

appears in Eq. (26) this is only a convenient way of writing the curl of an 
H field times a;o/(i£o), and variational error in E has no direct effect in the 
equation. 

If first order perturbation is applied to the H-ficld eigenvalue equation, Eq. (23), 
we also arrive at the result given by Eq. (22). We note that this is only appar- 
ently contradicting Eqs. (2) and (5) of Ref. [7], since these equations present 
the first order perturbation results in Ae and A(i), respectively. 



5 Example: Effective index of a photonic crystal fibre 



We consider the photonic crystal fibre first studied by Knight et al. [14] in 
which air-holes (of diameter d) are arranged in a micron scale triangular lattice 
(with pitch A) and the core is formed by a defect consisting of a "missing" air 
hole in the otherwise periodic structure. From the scale invariance of Maxwell 
equations [1] it follows that the dispersion properties are uniquely defined in 
terms of normalized quantities, such as the normalized free-space wavelength 
A/ A and the relative air-hole diameter d/h.. In this application we investigate 
how changing the size of the relative air-hole diameter changes the normalized 
frequency, a;A/(c27r) — A/A. 

One might speculate whether such problems can be addressed with the aid of 
first-order perturbation theory. However, for silica the variations of the air- 
hole diameter amounts to index variations of the order 5e ^ £si ~ ^Air — 
2.10 — 1.00 = 1.10 so that the relative change is certainly not close to zero as 
assumed in the derivation of Eq. (22). Thus this problem certainly calls for ap- 
proaches beyond ordinary perturbation theory. Here we apply the variational 
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approaches suggested above. 



5.1 Numerical model 

In fibre optics it is common to consider the frequency dependent effective 
index, 

, , c3(uj) 

n,s{uj) = ^^, (27) 

rather than the usual dispersion relation u{P). For obtaining a reference basis 
set we employ a plane-wave method in a super-cell configuration [16] , to pro- 
vide us with the cigcnmodes at a give frequency. Using these we evaluate the 
frequency shifts with both the variational method in E and H and diagonalize 
the operators to find the lowest cigcnfrcquency. For the PCF the dielectric 
function is isotropic so e(r) = e{r) ■ Isxs where e(r) is a real scalar. 

A PCF has a two fold degenerate fundamental mode corresponding to the two 
polarization states in free-space and cylinder-symmetrical problems [17]. In 
the following we will limit ourself to a Hilbert space consisting only of these 
two modes {N — 2). The PCF considered is made of a triangular structure 
of holes with the relative diameter d/A ~ 0.48, see Fig. 1. We consider a 
normalized propagation constant of j3A = 41.4, and a refractive index for 
the silica of nsi = 1.45. These structure parameters correspond to 1550 nm 
wavelength operation for a typical large-mode area photonic crystal fiber [18] 
and should be considered fairly standard for PCFs. 

5.2 Numerical results 

For the numerical implementation we employ a cell of 600 x 600 x 1 = 36000 
grid points. The minimal length between the cores in the super-cell images 
is 6 A. For comparison we do a full calculation for each hole size. In these 
calculations we set the propagation constant /3 in Eq. (27) to a fixed value 
and calculate the corresponding frequency. 

The calculation time for the full calculation with the plane- wave code is 25x 10 
minutes (6 hours) on a dedicated 1.4 GHZ machine. The calculation time for 
the presented method is 10 minutes for a single full calculation, 3 minutes for 
outputting the dielectric functions with the same program, and 1 minute for 
calculating the frequency shifts. The program was run on a desktop 866 MHZ 
machine and file I/O consumed the larger part of the time. The dielectric 
functions used in the presented method were the same as were used in the full 
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Fig. 1. Cross section of a photonic crystal fiber with a relative air-hole diameter 
d/A = 0.48. Light can be localized to a core region, which is formed by leaving out 
a single air hole in the otherwise periodic array of air-holes. The contours show the 
real part of a transverse component of the E-field amplitude calculated with the 
plane- wave method for /3A = 41.4. 

calculation, to give maximal consistency between the results of the variational 
method and the full calculations. 

In Fig. 2 the calculated effective index versus normalized air-hole diameter is 
presented. It is seen that the error is of second order in the hole size. This is 
due to the variational principle Eqs. (15,26) which causes the frequency u to 
be overestimated leading to an ries ~ uj~^ that is too small. 



5.3 Discussion 

The cause of the error in our variational methods is the discrepancy between 
the exact eigenmodes and the eigenmodes of the reference system. 

The electrical and magnetic fields are rapidly decaying inside the air holes with 
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Fig. 2. Effective index versus normalized air-hole diameter of the fundamental mode 
in a photonic crystal fiber calculated for /3A = 41.4. The upper curve shows results of 
a plane-wave simulation while the three other curves show the presented variational 
approaches taking the plane- wave simulation for d/A = 0.48 as a starting point. 

an exponential tail. So when the air hole boundary is moved a considerable 
distance, the reference fields deviate considerably from the exact fields. The 
condition for an electric charge free system, V • [£(r)E(r)] = 0, is not fulfilled 
for the 'new' system in the E variational method, whereas the H field is not 
subjected to the same condition in the H variational method, since we have 
V ■ [/iH(r)] = . This is probably the cause to the H variational method 
being slightly closer to the full calculation in Fig. 2. The reference fields in 
both variational methods thus do not describe the actual physical problem, but 
since we are interested only in eigenfrequencies and not their corresponding 
fields, this docs not pose a problem since the error in eigenfrequencies is of 
second order in the field. 

One might speculate why first order perturbation theory seems to perform bet- 
ter than the variational method as indicated by Fig. (2). This can be explained 
by observing Eqs. (21) and (22). In our problem the numerically evaluation 
of the matrix element in the first line Eq. (21) is an overlap integral (over 
the area given by e(r) — eo(r) 7^ 0) times a prefactor which is proportional to 
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£:o(r)/(£o(r) + Se{r) — 1 which is constant in the e{r) — eo{r) ^ area. . Since 
integrals for first order perturbation theory and the variational approaches 
evaluate to the same value, but the prefactor is — ^^(r) in the second line of 
Eq. (21), leading to first order perturbation theory to spuriously correct the 
variational error, by over- and underestimating the frequency for dj K < 0.48 
and d/A > 0.48, respectively. Theoretically, first order perturbation theory is 
not correct to any order for high-index contrast systems. 

One could include a couple of extra eigenmodes for an improved basis, but be- 
yond the two fundamental modes, the following modes poorly describe changes 
due to varying hole size. We could have included all 600 x 600 x 1 = 360000 
eigenmodes in the variational calculation, but this would be equivalent to a 
full calculation, and thus would make only little sense. One could, though, 
include eigenmodes that are confined around the hole boundaries for an im- 
proved description. This would give a much better estimation of the frequency 
changes. 

Finally, it is also possible to reduce the calculation time of the reference system 
by using a coarser grid, and then afterwards interpolate the fields to a higher 
resolution. 



6 Conclusion 



We have presented a method that can calculate frequency shifts for arbitrary 
changes in the refractive index, both isotropic and anisotropic, without com- 
plications due to degeneracy of the eigenfrequencies. We have demonstrated 
that the method can estimate frequency shifts which are exact to first order 
in the basis set. 

The method is capable of calculating frequency shifts which are too small to 
be calculate efficiently with a full calculation. The theory of the method is 
transparent and provides a sound basis for more elaborate numerical models. 

The fist of potential applications of the presented theory includes fiber grat- 
ings, birefringence of the doubly degenerate fundamental mode in PCFs due 
to structural inhomogeneities [19], and speed up for structural topological 
optimizations of electromagnetic systems [20]. 
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